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1. Introduction 

We study the non-linear diffusion of a population of Cosmic Rays (CRs) after they escape from 
the acceleration site (for example, a supernova remnant shock). As a zero order approximation the 
problem can be described by an impulsive injection of CRs of a given energy at a given location 
in the Galaxy. We study the transport of CRs taking into account the growth of Alfven waves 
caused by particle-wave interactions, and their back-reaction on the CR diffusion. We also consider 
possible linear wave damping mechanisms that limit the growth of the waves. The problem is 
described by two differential equations that govern the evolution of the CR pressure and of the 
wave energy density as a function of time and distance from the acceleration site. We present 
numerical solutions of these two coupled equations for different values of the physical parameters 
involved in the problem. 


2. The model 


In our model the transport of CRs is assumed to be regulated by the resonant scattering off 
Alfven waves, i.e. a CR of energy E resonates with waves of wave number k= \/ri{E) where 
is the particle Larmor radius (see e.g. ||T]]). The (normalized) energy density I{k) of Alfven waves 
is defined as: 

^ = |//(i)dl„i, (2.1) 

where So is the ambient magnetic field and 8B the amplitude of the magnetic field fluctuations. 

According to quasi-linear theory, CRs diffuse along the magnetic field lines with a diffusion 
coefficient equal to |0 : 


4 c rt{E) _ Db{E) 
3711{k) l{k) 


( 2 . 2 ) 


which can be expressed as the ratio between the Bohm diffusion coefficient Db{E) and the energy 
density of resonant waves l{k = \/ri). Formally, quasi-linear theory is valid for 8B/Bq <C 1. In 
this limit the diffusion perpendicular to the field lines can be neglected, being suppressed by a 
factor of {8Bk/B()y = I{k)^ with respect to the one parallel to Bq (see e.g. [^,@1). This implies that 
under the condition that 8B/Bq remains small, the problem is one-dimensional. 

We consider the streaming of CR to be the main source of Alfvenic turbulence. For CRs that 
stream along the direction of Sq (that we consider aligned along the coordinate z) the growth rate 
of Alfven waves Fes is proportional to the product between the Alfven speed Va and the gradient 
of the pressure of resonant CRs, and can be defined as (see e.g. iJ): 


3-TcrI 


-Va 


dPcR 

dz 


(2.3) 


The sign indicates that only waves traveling in the direction of the streaming are excited. For 
dimensional consistency, the pressure of CRs Pcr has to be normalized to the magnetic energy 
density Bq/Sti. 

On scales smaller than the magnetic field coherence length, a flux tube characterized by a 
constant magnetic field strength So and directed along the z-axis can be considered. Two coupled 
equations (for the evolution of CRs and waves along the flux tube) must be solved: the diffusion 
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coefficient of CRs of energy E (Eq. 2.2), indeed, depends on the energy density of resonant waves 


I{k), and in turn the growth rate of these waves (Eq. |2.3|) depends on the gradient of the pressure of 
resonant CRs. The two coupled equations then read: 

^ f Db dPcR 
dz 


dPcR , dPcR 
■ + Va- 


dt 


dz 


I dz 


(2.4) 


dl dl dPcR 


(2.5) 


where the left side in both expressions is the time derivative computed along the characteristic of 
excited waves: 

d ^ 

( 2 . 6 ) 


d _ d d 
dt dt~'^ ^ dz 


The advective terms VAdPcR/dz and VAdljdz are neglected in the following since we verified that 


they play little role in the situation under examination. The last two terms in Eq. 2.5 account for 
possible mechanisms of wave damping, operating at a rate E^, and for the injection Q of turbulence 
from an external source (i.e. other than CR streaming). In this work we consider only linear 
damping mechanisms, for which the damping rate is independent of space and time. The term 
representing the external source of turbulence can be set to Q = ir^Io, so that when streaming 
instability is not relevant the level of the background turbulence is at a constant level I = /q. 

As already pointed out by [|^, the approach described above decouples the process of acceler¬ 
ation of particles, which operates, for example, at a SNR shock, from the particle escape from the 
acceleration site. Though such a separation might seem artificial, the problem defined above can 
still be useful to describe the escape of particles form a dead accelerator, in which the acceleration 
mechanism does not operate any more, or operates at a much reduced efficiency |^. This situation 
would probably apply to the case of old SNRs. On the other hand, suggested that Equations. 2.4 


and 2.5 could be also used to describe an intermediate phase of CR propagation in which the CRs 


have left the source but are not yet completely mixed with the CR background. Eor the case of a 
SNR shock, the equations above could thus be applied to those CRs characterized by a diffusion 
length large enough to decouple them from the shock region. Typically, this happens during the 
Sedov phase to the highest energy particles confined at the SNR shock when the diffusion length 
Dr/U s gradually increases with time up to a value larger than some fraction % of the SNR shock 
radius Rs, where Ug is the shock speed and x ~ 0.05...0.1 


In both scenarios, the initial conditions for Eqs 2.4 and 2.5 can be set as follows: 


PCR = PcR 
PcR = 0 


and 

and 


iKil 

/ = /o«l 


for z<a 

for z> a 


(2.7) 

( 2 . 8 ) 


where a is the spatial scale of the region filled by CRs at the time of their escape from the source. 
Eollowing 0, we introduce the quantity IT, defined as: 


n = 


Db 




CR 


where: 


= 


noo 

/ dz PcR = PcrU 
Jo 


(2.9) 


( 2 . 10 ) 
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To understand the physical meaning of IT, consider the initial setup of the problem, in which CRs 
are localized in a region of size a. The CR pressure within a is and then ^CR = Pcr^- The 
initial diffusion coefficient outside the region of size a is equal to Db/Iq- The problem of transport is 
characterised by three relevant timescales: i) the growth time Tg ^ {VA/IodPcR/dz)^^ ~ alQ/VAP^R, 
ii) the time it takes the CR cloud to spread due to diffusion Xdiff ~ jD o^Iq/Db, and iii) the 
damping time = 1 /2Td. In order to have a significant growth of waves due to CR streaming, the 
timescale for growth must be shorter than the two other timescales: Xg < min(T^,yy, T^). In terms of 
the parameter IT this conditions reads: IT > max(I, Zdiff /id). It is evident, then, that the parameter 
n regulates the effectiveness of CR-induced growth of waves. In the absence of a damping term for 
waves, n > I is a necessary condition for streaming instability to be relevant, while in the presence 
of efficient wave damping, a more stringent condition on IT applies. 

For n < max(l,Td( 7 //Td) CRs play no role in the generation of Alfven waves, and Equa¬ 


tion 2.4 can be solved analytically: 


PcR = 


h 


TIDbI 


fi>CR e 




( 2 . 11 ) 


Equation [2.1 1| is referred to as test-particle solution. When wave growth cannot be neglected, 
the solution deviates from this analytic test-particle solution. In the extreme scenario IT >> 


max(l,Td 7 y/Td) waves grow so quickly that the the diffusive term in Eq. ^ becomes negligi¬ 
ble when compared to advection term VAdPcR/dz, and the advection term cannot be neglected 
anymore. This describes a situation in which CRs are "locked" to waves and move with them at a 
velocity equal IoVa [W]. An identical result was found by [|TI|] in a study of the escape of Ri MeV 
CRs from sources, and also suggested by [|T^]. 


3. The method 


In order to solve Eqns. 2.4 and 2.5 it is convenient to perform the following change of coordi¬ 


nates 


s = 


T = 


and re-normalise the parameters as follow: 




r' = 


In these notations (and after dropping the advection terms) the equations become: 

d^CR d ( 1 B^cr' 


dz 

dW 

dx 


ds 


W 
dS^CR 
ds 


ds 

-2r'(w-Wo) 


(3.1) 

(3.2) 

(3.3) 

(3.4) 


Their solution (i.e., ^ and IT as a function of the variables T and s) depends only on three 
parameters: the initial values and Eq. Note that: i) ^CR ~ ii) ^0 = VaoIDq, and iii) 

r' = Eq (since we consider here only a linear damping constant in both time and space). In these 
notations, the condition for effective growth of waves is IT > max(l,r'ITo). 
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Figure 1: Left: ion-neutral damping rate (L^^, solid curves). Farmer & Goldreich damping rate (Fj*^, dashed 
curves) and ratio between Alfven speed and initial size of the source iVAja, dot-dashed curves), as a function 
of the CR energy £. Different colours refer to two different ISM phases: warm ionised medium (WIM, red), 
and warm neutral medium (WNM, blue). Right: parameter If as a function of the CR energy for a WIM 
and WNM (solid lines). In absence of efficient wave damping, values of 11 < 1 correspond to a test particle 
solutions, since the wave growth rate is inefficient. Conversely, for 11 larger that ^10 the growth rate 
is large and the diffusion coefficient becomes smaller than the Bohm diffusion coefficient Db- In presence 
of efficient wave damping, both these limits are shifted to higher values, depending on the value of F^. 


We present examples of numerical solutions to the coupled equations and for values 
of the parameters in the ranges: IT = [2 — 20], Wb = [10^^ — 10^^], and F' = [0 — 100]. In order 
to understand why we chose these values, and to which physical conditions they correspond, it is 
necessary to specify the values of VA,a,Do, Db, and the damping mechanisms. These quantities in 
general depend on the properties of the ISM, and on the particle energy. We consider two different 
ISM phases: a warm neutral medium (WNM) and a warm ionised medium (WIM). The typical 
values of density, temperature, ionisation fraction, and magnetic field strength for these two phases 
are taken from [13]. In these ISM phases, the two most relevant linear damping mechanisms are the 
ion-neutral friction (F^^) and the damping by background MHD turbulence suggested by Farmer & 
Goldreich [ 14] (F^^). The value of F^^ as a function of the particle energy is shown in Fig. (left 
panel, long-dashed lines), for both ISM phases. The ion-neutral damping rate is estimated by nu¬ 
merically solving the equations in [^]. The results are shown in Fig. |I] (left panel, solid lines), for 
both ISM phases: at low energies F^^ = V/a?/ 2 (where V/at is the ion-neutral collision frequency), 
while at high energies F^^ oc The level of coupling between ions and neutrals affects also 
the Alfven speed, that is in the range 10^ — 10^ cm/s. Since it is useful for our calculations, we 
show in Fig. |I]the value of the ratio Va/a, as a function of energy (dash-dotted lines). In order to 
estimate a we consider a « Resc{E), where Resc{E) is the escape radius and depends on the CR 
energy. We model its dependence on the energy following [1^]: we found that its value ranges 
between a few pc (for the highest energy particles) and ~20 pc (for ~GeV particles). Since we are 
interested in following the CR diffusion from the scale a where they are released up to a distance 
z < 10^ pc, this roughly corresponds to solve the normalised equations in a range of s from 1 to 
several hundreds. Moreover, we want to study the CR pressure at different times t in the range 
10^^ yr to < 10^ yr. Considering the value of Va/a as a function of energy displayed in Fig. |^, this 
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corresponds toT = 10 ^ — 1. 

In order to identify the relevant values of the parameters IT, Wb,r' we proceed as follow: 

• r'- the damping rate Fj = mar(r^^,r™) ranges between 10^^ and 10^^^s^^ (Fig. |I], left 
panel), depending on the particle energy and ISM phase. This corresponds to i<r<io4. 

• n - the values of IT as a function of the CR energy, for the two different ISM is shown in 
Fig. [T] (right panel, solid lines). 

• Wb - for the background diffusion Dq we can refer to the average Galactic value Dq = Djsm = 
10^^(£'/10GeV)°-^cm^/s. Since Fra ranges between lO^^cm^/s and lO^^cm^/s (at low and 
high energies, respectively), the range of relevant values for Wq is 10~^ < Wb < 10~'. 


4. Results and Discussion 


We numerically solve Eqs. and and show the results for £P[s, t) vs. s (solid lines in 
Fig. 1^ at three different normalised times T = 10^^, 10^^, 1 (corresponding to different colours in 
the figure). The test particle solution is also shown for comparison (dotted lines). Each row in the 
figure corresponds fo resulfs obfained by varying only one parameter af fhe lime. 

The upper row shows fhe impacf of fhe paramefer IT for fixed values of F' = 10 and Wq = 10^“^. 
Since for IT < 1 fhe numerical solulion coincides wilh fhe lesf paiticle solulion, we show our resulfs 
for values of IT > 1. When IT = 2, fhe solulion slarls lo deviale from fhe lesf parlicle case: due lo 
the growth of wave by streaming instability the diffusion is slower, the CR pressure in the vicinity 
of the accelerator is larger, and the diffusion length is smaller. The effect is more and more evident 
for increasing values of IT (upper row, middle and right panels). 

The middle row shows the effect of wave damping, for fixed IT = 10 and Wq = 10^^. Since 
n >> 1, when wave damping is negligible (F'=0, lefl panel) fhe diffusion is suppressed as com¬ 
pared lo fhe lesf paiticle case. Earger values of F' limil fhe growlh of Alfven waves and fhe solulion 
approaches fhe lesf parlicle solulion. Note lhal F'=10 (middle row, cenlral panel) corresponds lo a 
damping mechanism lhal becomes imporlanl al T > 0.1. For Ihis reason Ihe solulion al T = 0.01 
and T = 0.1 are unaffected by Ihe increased damping rate, and only Ihe solulion al late lime T = 1 
is modified, and approaches Ihe lesl parlicle solulion. The righl panel shows Ihe case of a even 
faster damping mechanism, F'=100, whose effecl slarls lo be relevanl al limes corresponding lo 
T > 0.01, explaining why Ihe solulion al intermediate and late limes approaches Ihe lesl parlicle 
solulion, while Ihe solulion al early lime is unaffected by wave damping. 

The bottom row shows calculalions performed al IT = 5 and F'=10, wilh Wb varying from 10~^ 
lo 10^^. A large Wq (lefl panel) corresponds lo slow diffusion and small diffusion lenglh, both in the 
test particle case and in the numerical solution. The two solutions differ at early and intermediate 
times because the growth rate is faster than the diffusion rate, and suppresses the diffusion. At later 
times instead the growth rate decreases and the numerical solution does not appreciably differ from 
the test particle case. A change in the initial value of the background turbulence has the same effect 
on both the analytical and numerical solution: the diffusion is faster for decreasing values of Wq 
(see bottom row, from left to right), and the CRs can reach larger distances. The deviation from a 
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n=2 n =6 n=20 



s s s 

r'=o r'=io r'=ioo 



s s s 

Wo=10" Wo=10'^ Wo=10'^ 



s s s 


Figure 2: Normalized CR pressure as a function of the variable s, which represents the distance from the 
source z normalised to the size of the source a. In each row the results are shown by varying only one 
of the three parameters (whose value is reported above each panel) and keeping the other two fixed to the 
values reported in the shaded area. The solid curves show the results from the numerical computation, while 
dotted curves represent the test particle solution, where the role of streaming instabilities in enhancing the 
turbulence is neglected. Different colours refer to different normalised times T = 10^^, 10^ \ 1. 
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test particle solution however remains unchanged, since both the growth rate and the diffusion rate 
scale as Wq. 

This study outlines the importance of considering, in the physics of CR transport, both the 
growth of waves due to streaming instability and the role of mechanisms that damp the turbulence. 
When damping is neglected, even relatively small values of the integrated CR pressure IT = 10 
(corresponding, for example, to particles of a few TeV, in a WIM - see Fig. [I]) in a environment with 
small background turbulence Wq = 10^"^ (relevant for TeV particles when the background diffusion 
coefficient Dq is similar to the average Galactic diffusion coefficient) lead to the conclusion that 
diffusion is strongly suppressed by the efficient growth of waves, and particles cannot travel large 
distances, even at late times. This situation is shown in Fig. (middle row, left panel) where the 
numerical computation is performed in absence of a wave damping term: the diffusion is strongly 
suppressed also at T = 1, corresponding to old remnants (tage rsj 10^ yr). 
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